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^ Abstract 

o 

We describe a method for detecting, locating and quantifying sources of gas emissions to the at- 
^ i mosphere using remotely obtained gas concentration data; the method is applicable to gases of 

environmental concern. We demonstrate its performance using methane data collected from air- 
craft. Atmospheric point concentration measurements are modelled as the sum of a spatially and 
7—i temporally smooth atmospheric background concentration, augmented by concentrations due to 

local sources. We model source emission rates with a Gaussian mixture model and use a Markov 
random field to represent the atmospheric background concentration component of the measure- 
ments. A Gaussian plume atmospheric eddy dispersion model represents gas dispersion between 
sources and measurement locations. Initial point estimates of background concentrations and source 
emission rates are obtained using mixed ^2-^1 optimisation over a discretised grid of potential source 
locations. Subsequent reversible jump Markov chain Monte Carlo inference provides estimated val- 
ues and uncertainties for the number, emission rates and locations of sources unconstrained by a 
grid. Source area, atmospheric background concentrations and other model parameters are also 
estimated. We investigate the performance of the approach first using a synthetic problem, then 
apply the method to real data collected from an aircraft flying over: a 1600km^ area containing 
two landfills, then a 225km^ area containing a gas flare stack. 
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1. Introduction 



There is growing interest in developing methods for detecting and locating sources of gas emissions 



into the atmosphere. Greenhouse gases are of intense interest (e.g. Chen and Prinn (2006); Shak- 



hova et al. (2010)). Other applications include monitoring toxic gas emissions, locating explosives 



from their volatile emissions (e.g. Bhattacharjee ( 2008| )]), mapping naturally occurring gas seeps 
for oil and gas exploration (e.g. Hirst et al. (2004)), identifying sources of nuisance odours, and 
even understanding how moths are able to find mates by detecting pheromones at concentrations 
corresponding to individual molecules (e.g. Vergassola et al. (2007)). For greenhouse gases and 



oil and gas exploration the goal is to locate sources and quantify emission rates. For explosives, 
nuisance odours and moths, locating the source is sufficient! 

In this paper we concentrate on the task of detecting, locating and quantifying the emission rates of 
sources of a single gas species of interest. While the method we have developed is broadly applicable 
to any gas or passively transported, detectable atmospheric component (e.g. aerosols, radon, smoke, 
dust, viruses) we concentrate on methane emissions using data acquired during development and 
testing of an airborne system for mapping natural gas seeps for use in hydrocarbon exploration over 
areas typically up to 5000 km^ per flight. We apply the method to measurements from test flights 
around two modern Canadian landfills and a flare stack within a modern natural gas processing 
facility in North Africa, collected at ranges of up to 12km downwind of sources. We make inferences 
about source emission rates using gas concentration measurements. We measure concentrations by 
volume and use our understanding of gas dispersion to relate these to source emission rates expressed 
as m^s~^ of pure gas per source. For area sources emission rates can also be expressed as mass flux; 
i.e. mass emission rate per unit time per unit area. Critical to achieving this goal is estimating the 
level of atmospheric background concentration so we can identify the additional concentration over 
and above background, attributable to the local source(s) of interest. The novelty of the current 
work lies in the tailored application of a combination of standard statistical modelling components 
and inference tools for inversion in remote sensing. 



The extensive literature on inversion and the related field of compressive sensing (e.g. Donoho 



(2006)) includes contributions from the statistics, applied mathematics, electrical engineering and 



physics communities. Sambridge and Mosegaard (2002) reviews the development and application 



of Monte Carlo methods for inverse problems in the Earth sciences. Rao (2007) reviews source es- 



timation methods for atmospheric releases of toxic agents, including forward modelling (potentially 
using Bayesian inference) and backward transport modelling, emphasising the need to assess uncer- 



tainties in characterisation of sources using atmospheric transport and dispersion models. Senocak 



et al. (2008) discusses early detection of the location and size of a contaminant release into the 



atmosphere from a network of environmental sensors using a Gaussian plume forward model with 
stochastic parameters and Bayesian inference using Markov chain Monte Carlo (MCMC). It is noted 



that the main distinguishing feature of Bayesian inference as opposed to optimisation (e.g. Thom- 



son et al. (2007)) for source estimation is that the former estimates probability distributions for 
parameters of interest and quantifies the uncertainty in the estimated parameter, whereas the latter 
provides point estimates for the parameters of interest through optimising an objective function. 



Keats et al. (2007) note that determining the source of an emission from the limited information 



provided by a finite and noisy set of concentration measurements obtained from real-time sensors 
is an ill-posed inverse problem. They show that solving the adjoint advection-diffusion equation 
just once per detector location allows efficient forward model estimation and Bayesian inference 



using MCMC. Long et al. (2010) uses a genetic algorithm to find the combination of source loc- 



ation, height and strength, surface wind direction and speed, and time of release that produces a 
concentration field that best matches sensor observations. A rationale is developed to specify the 
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minimum number of sensors necessary to estimate the source term and to obtain the relevant wind 



information to a given precision. Rudd et al. (2012) presents an inverse modeUing technique to 



estimate source strength and location, together with the uncertainty in those estimates, using a lim- 
ited number of measurements from a sensor network. Experimental design aspects are addressed, 
including the optimal number and configuration of sensors for a given measurement campaign, and 
the minimum period of observation for source detection with confidence. The need to relate un- 
certainty in estimated source properties to those of the input data is emphasised. [Gyarmati-Szabo 



et al. (2011 ) introduces a Bayesian multiple change-point model for monitoring of air quality stand- 
ards by pollutants such as nitrogen oxides and carbon monoxide. Change-points are identified and 
the rate of occurrence of air quality threshold exceedence estimated using a reversible-jump MCMC 
approach. 

The layout of the paper is as follows. Motivating landfill and flare stack applications are first 
described in section [2} Section [3] outlines the modelling procedure and illustrates it in application 
to a synthetic problem. Application of the method to the landfill and flare stack measurements 
is described in section |4j Section [5] provides a discussion of findings and suggestions for future 
development. Modelling details are relegated to three appendices, describing background modelling 
(|A|), initial parameter estimation using mixed optimisation ([B|) and mixture modelling ([C|. 



2. Data 



We used an ultra-sensitive, high precision methane gas sensor, mounted in an aircraft, to measure 
a continuous stream of air from the leading edge of a wing - well away from any fuel/lubricant 
or engine exhaust fumes. The sensor continuously measures gas concentration and passes data 
to the aircraft's data logging system together with GPS, radar altitude, barometric pressure, air 
temperature and wind velocity data, and several system control parameter values. The sensor 
delivers better than 1 ppb (part per billion by volume) precision concentration data with a response 
time of approximately 1 second. Flight data are subsequently merged with specialist meteorological 
data, including additional wind, atmospheric boundary layer depth and auxiliary data: such as the 
air sample transit time from sample inlet to sensor measurement chamber. 

The data sets presented here are atypical of surveys aimed at hydrocarbon exploration, in that 
we know the source locations within in the survey areas. Consequently, these data sets provide a 
valuable test of our measurement and analysis procedures. Source related concentrations are up 
to two orders of magnitude greater than those typically encountered in natural seepage surveys. 
Landfill measurements provide a direct test of the total system performance; flare stack measure- 
ments provide a means of determining the gas sample transit time in flight at operational speeds; 
this is essential to correctly assign concentration data to measurement locations. 



2.1. Landfills 



Atmospheric methane is responsible for about one third of the global warming effect of CO2 des- 
pite CO2 concentrations being approximately 220 times greater than those of methane. There 
are strong incentives to reduce methane emissions to the atmosphere. Landfills are prodigious 
sources of anthropogenic methane; about 25% of United States' methane emissions are from land- 



fills, (International Energy Agency, 2009). 



The global average atmospheric methane concentration is approximately l,820ppb, increasing at 



approximately 8ppb per annum (Dlugokencky et al. , 2009). Local atmospheric background methane 



concentration can vary by approximately 20 ppb during daytime due to changes in Atmospheric 
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Boundary Layer (ABL) depth. This layer of the atmosphere effectively contains all ground level 
emissions for that day. Its growth is driven by solar heating of the ground. In effect, ground level 
emissions accumulate near the ground during the night when the ABL is thin and are diluted into the 
growing volume of the ABL during daytime. It is important to account for the associated changes 
in local background concentration, so as to more precisely determine that portion of measured 
concentration attributable to the local sources of interest: since source related concentrations can 
be comparable to the much longer term changes in background. 

The top of the ABL constitutes a "ceiling" on vertical transport of gases from the ground, and 
reduces the rate of dilution with downwind distance from the source. This effect must be included 
in the gas dispersion model used to relate measured concentration to source strength. It can be 
modelled as the sum of multiple reflections from the ABL "ceiling" and ground surface. At longer 
ranges such as those prevailing here, there will be multiple bounce terms successively from the ABL 



"ceiling" and the ground surface, until the air within the ABL is well mixed (Gilford (1976)). The 



aircraft must be within the ABL if it is to detect concentrations from ground level sources. 

The data presented here correspond to measurements at the aircraft flies at approximately 200m 
AGL (above ground level) at a speed of approximately 50ms~^; ABL depth is approximately 400m. 
The area of interest is 40km x 40km and the flight time approximately 80 minutes. Figure [T] shows 
the flight track of the aircraft around each site, tracing a loose serpentine pattern downwind of each 
landfill. Initial average wind speed and direction (based on measurements near the site and model 



data provided for multiple altitudes by the UK Meteorological Office, UKMO Davies et al. (2005) 
) are about 6.5ms^^ and 033° degrees meteorological. (Wind direction is defined as the direction 
from which the wind blows, in a clockwise sense, with North as 0°. In this case, the wind direction 
is approximately North-Easterly, i.e. air moves approximately to the South West.). 

Figure [T] shows the flight trajectory around and in the vicinity of the two landfills; Figure [2] shows 
measured methane concentrations in time. Inspection of the concentration trace along the tra- 
jectory in Figure [T] indicates the wind direction is approximately constant for the period of the 
flight. 

2.2. Flare stack 

The flare stack flight comprises 8 separate multi-looped circuits of the flare at altitudes from 150m 
to 350m AGL; ABL depth is greater than 1500m. Each circuit intersects the plume three times 
at different angles and ranges to probe the 3-D structure of the dispersion plume. Additionally 
measurements are used to determine the air sample transit time to provide consistent plume de- 
lineation. The flare stack is 50m high, situated within a recently completed natural gas processing 
plant at a coastal location, where winds are variable, probably burning methane and light hydro- 
carbons. Photographs of the flare stack show it to be a clean yellow flame, which suggests it is 



burning at high efficiency (Kearns et al. , 2000). Methane content of a few percent is expected in 



the residual unburnt gas; this is due to thermal decomposition and incomplete combustion of the 



fuel- rich hydrocarbon feed ( Pekalski et al. , 2005 ) . 



Initial UKMO model-based average wind speed and direction are llms^^ at —243°. Figure [s] shows 
the flight track around and in the vicinity of the flare stack; Figure [4] shows measured methane 
concentrations against time. Inspection of the concentration trace along the track in Figure [3] 
suggests that the UKMO predicted wind direction is inaccurate for the portion of the flight: which 
is late in the afternoon and situated over the coastline. 
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Figure 1: Flight track around and in vicinities of two landfills. Blue marker size and colour saturation indicate strength 
and location of measured methane concentrations. Arrow shows average direction of predicted air movement during 
flight. Polygons show perimeters of methane-emitting landfill areas. Dimensions in km. Aircraft takes off in North 
East corner and flies over Westerly landfill first. 
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Figure 2: Methane concentrations along landfill flight path as a function of time. 
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Figure 3: Flight track around and in vicinity of the flare stack. Blue marker size and colour saturation indicate strength 
and location of measured methane concentrations. Arrow shows average direction of predicted air movement during 
flight. Black annulus indicates location of flare. The visible discrepancy in alignment of significant concentrations 
and flare stack with respect to arrow indicates error in predicted direction of air movement. Aircraft enters from the 
SW corner and leaves NE corner. 

3. Model 

3.1. Model form 

We seek to locate and quantify emission rates of methane sources given n observed atmospheric 
concentration measurements y = along airborne trajectory x = {xj}"^]^. The overwhelming 

majority of concentration measured is the atmospheric background contribution of about ISOOppb 
with an additional — lOOppb (typically) attributable to local ground level sources. Background 
level varies in space and time. Poor estimation of background concentration disproportionately im- 
pacts estimation of concentration attributable to local sources. Hence we prefer to jointly estimate 
background and source contributions, y is modelled as the sum of a slowly varying background 
b = along the trajectory and contributions due to a distributed group of m sources at 

ground level locations z = {zj}JL^ with emission rates s = {sj}JLi and auxiliary characterist- 
ics C = {Cj}^^^. Measurements along the trajectory are assumed to be made with independent 
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Figure 4: Methane concentrations along flare stack flight path as a function of time. 



identically-distributed additive Gaussian errors e = {ei}^^^, with ~ N{0,a'^). Steady-state gas 
transport between a source of unit emission rate at location Zj and measurement location Xj in 
wind field W is given by coupling function aij = a(xj, Zj, Cj; W). With A = we adopt 

the model: 

y = As + b + e (1) 

Model mis-specification can be diagnosed by analysis of residuals, e.g. the elements of e in equation 
[T] should represent a random sample from a normal distribution with mean zero and constant 
variance. 



Plume model 

The wind field W is described by wind vector U at measurement location x, and horizontal and ver- 
tical plume opening angles and 7y respectively. We approximate coupling a{x, z,w;\J ,^Hi7v) 
between a unit source of half width w at location z and measurement location x in wind field W 
by a Gaussian plume model. The measurement location relative to the plume is expressed in terms 
of downwind distance 5r, and horizontal and vertical offsets 5h and 5v of measurement location 
with respect to wind vector, an = 5_Rtan(7//) + w and ay = (5/jtan(7y) play the role of plume 
"standard deviations" in horizontal and vertical directions respectively. 

where D is the maximum altitude of the atmospheric boundary layer above ground level (referred 
to as the atmospheric boundary layer depth) . The sum of four exponential terms represents plume 
reflections in the ground and at the interface between atmospheric boundary layer and free atmo- 
sphere above it at altitude D. Values of U, D, and 7y are obtained directly from wind field 
data supplied by UKMO. The plume model parameters are illustrated in Figure [5] Typical plume 
characteristics are shown and discussed in Section [4.11 
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Figure 5: Illustration of plume model parameters. Red line: source height H . Magenta line: downwind distance 5r 
of measurement location relative to the source. Cyan line: horizontal offset 5h- Green line: vertical offset Sv- Blue 
line: source half width w. Thick black horizontal lines perpendicular to wind direction represent ground level and 
top of ABL at height D. The current plume model allows for reflections from ground and ABL ceiling. In this figure 
Sn is too short for reflections to be effective. The hatched blue area represents the marginal variation with 5h of the 
value of a for Sv =Q and fixed Sr. The hatched red area represents the marginal variation with Sv of the value of a 
for fixed 5h and Sr. The locus of the black dot (corresponding to a contour of constant a for fixed Sr) is drawn as a 
black closed curve. 

Background model 

Well-mixed background gas concentration b along the trajectory is assumed to be positive and 
smoothly-varying spatially and temporally. We assume it can be represented by an appropriate 
set of r-dimensional smooth spatio-temporal basis functions: 

b = P/3 (3) 

for an n X r matrix P of bases, where /3 = {likYk=i parameters to be estimated. As a default 
approach, we assume P to be the identity matrix and adopt a Gaussian Markov random field model 
for /3(= b). In this case, /3 is a random vector with prior probability density function: 

/(/3) (X exp{-f (/3 - /3o)^(J/3)(/3 - /3o)} (4) 

with pre-specified precision matrix and tuning parameter /i. Due to the random field's condi- 
tional independence structure, J/3 is guaranteed to be sparse, allowing efficient parameter estim- 
ation. As outlined in Appendix we specify the precision matrix such that pairs of locations 



8 



aligned with tlie wind direction have higher dependence. We also consider other representations of 
the form b = P/3, e.g. polynomial or spline bases for which r <^ n (see Appendix [A|) . 



3.2. Initial parameter estimation 

Given measured concentrations y on trajectory x and associated wind field data ({U(xj)}[L]^, 7h, 7y), 
initial point estimates for source emission rates s and background parameters (3 are obtained as 
follows. We assume a Laplace prior for s with pre-specified precision Q and tuning parameter A: 

/(s)ocexp{-A||Qs||i} (5) 

Since the likelihood corresponding to model [T] is: 

/(y|s,/3) ocexp{-2^ps + P/3-yf}, (6) 

by applying Bayes theorem, the posterior parameter density becomes: 

/(s,/3|y)«/(y|s,/3)/(s)/(/3) (7) 

In particular the maximum a-posteriori parameters are obtained by maximising /(s,/3|y) with 
respect to s and /3, or equivalently by minimising — logg /(s, /3|y). In optimisation terms, initial 
parameter estimation can be stated as: 

argmin,^^ _l^Ps + P/3 - yf + f (/3 - /^of J(/3 - /3o) + A||Qs||i (8) 

where terms in fj, and A can be viewed as regularisations that impose background smoothness and 
source sparsity respectively. We further choose to restrict the domain of source elements such that 
Sj E [0, Smax], and background bj G [0, + r] for tolerance r. Details are given in Appendix [Bj 



3.3. Mixture model 



Full parameter estimation is performed using a mixture modelling approach. We assume that each 
of m sources can be represented as a two-dimensional Gaussian kernel located at Zj with half width 
(corresponding to the standard deviation of the Gaussian) and source emission rate Sj. Using 



reversible jump Markov chain Monte Carlo (RJMCMC) simulation (Green (1995)), we treat m 



as a random variable, and estimate the joint distribution of m and all other model parameters. 
We can also make inferences about apparent bias and/or uncertainty in wind field parameters and 
measurement error. Bias-correction of wind direction proves to be important in some applications. 



RJMCMC for mixtures of univariate Gaussians was considered by Richardson and Green (1997) 



and extended to multivariate Gaussian mixture models by Zhang et al. (2004). 



Markov chain Monte Carlo (MCMC) is a simulation procedure for Bayesian inference, which exploits 
the fact that Markov chains have stationary distributions exhibiting reversible transitions. We 
seek to create a Markov chain with the posterior distribution /({z, w, s, /3}|y) as its stationary 
distribution using the Metropolis-Hastings algorithm (jMetropolis et al. ( 1953 ) and Hastings ( 1970 )). 
Extending the posterior density in ([7| to include source locations and half widths, and writing the 
resulting parameter set {z, w, s, (3} as for brevity, the expression for the posterior parameter 
density becomes: 

f{e\y) cc f{y\e)f{e) (9) 

Enumerating the constant of proportionality in ^ is generally computationally costly, and hence 
so is direct sampling from the posterior. Fortunately MCMC (e.g. using Metropolis-Hastings 
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acceptance sampling or similar) circumvents the need to evaluate the constant of proportionality. 
We proceed by judiciously partitioning the set of model parameters 6 (to exploit problem structure 
and ensure reasonable MCMC performance), so that dependent parameters appear in the same 



subset 6^ of parameters indexed by k (see, e.g. Gammermann and Lopez (2006)), and that different 
sampling techniques (such as Gibbs sampling, see below) can be exploited for different subsets. The 
conditional posterior distribution of parameter subset 0^ becomes: 



fie^\y,ej,)^fiy\e^,ej,)f{e,\ej,) 



(10) 



where 0^ represents the remaining model parameters. The Metropolis-Hastings algorithm then 
provides a basis for acceptance of a candidate {O'^jO^} (with 6'^ generated by sampling from a 
suitable proposal (7(0^) ^kI^k); such as a multivariate Gaussian centred at 6^) given the current 
Markov chain position {6^, 6^} with probability 0(6^, 0'^\6k): 



f{e',\y,e^)q{ele,\e^) 



(11) 



Since the conditional posterior appears in both the numerator and denominator of (11), the fact 



that the conditional posterior is only known in (10) up to a constant of proportionality is of no 



consequence. Starting from an arbitrary starting point, having allowed for a period of burn-in (to 
facilitate convergence (Gammermann and Lopez, 2006)), the sequence of points 9 = {z, w, s, (3} on 
the Markov chain will converge to a (dependent) sample from /(z, w, s, (3\y). In this way, for a fixed 
number m of sources, we can estimate the joint posterior distribution of the model parameters. 

RJMCMC allows sampling from distributions for which the number of sources m (and hence the 
total number of model parameters) is not fixed. As explained in Appendix [C| the Metropolis- 
Hastings algorithm can be extended to accommodate "birth" of a new source, "death" of an existing 
source, coalescence of neighbouring sources and source division. The Markov chain will therefore 
explore estimates of z, w and s of different dimensions together with /3. Inference can be extended 
to include estimation of quantities such as the measurement error standard deviation a^, bias of 
wind vector U and horizontal and vertical plume opening angles 'Jh and 7v. In the analysis 
reported in section |4| (see ([l])), and additive wind direction bias (used in ([2])) are included as 
model parameters in the Bayesian inference. 

When the conditional posterior distribution f (9 i^\y , Oj^) can be written in closed form, values of 
9k can be sampled directly given current values of 9^. This approach, known as Gibbs sampling 



(e.g. Geman and Geman (1984)) or sampling from full conditionals, avoids the need for acceptance 
sampling. In the current work, the conditional distribution of background parameters (3 is known 
in closed form and is amenable to sampling from full conditionals. When full conditionals are not 
available we use the Metropolis-Hastings algorithm. 



The initial optimisation solution (section 3.2) is sampled to give a suitable starting point to accel- 



erate MCMC convergence to the stationary distribution. 



4. Application 

4.1. Illustrative analysis 

We illustrate the approach using a synthetic problem. Gas plumes from 10 randomly located ground 
level methane sources are generated with emission rate of O.lm^s"^ propagating in a variable 
wind field generated using a random walk, with wind speed and direction of [6.3, 6.6]ms~^ at 
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between [218,222]° meteorological, and plume opening angles of between 'jh = [11.5,13.9]° and 
7y = [11.5, 13.9]°, within an area of 40km x 40km. We simulate methane concentrations on a 
flight path of Ihour 9mins sampled at 3s intervals yielding 1379 observations, and add a constant 
background of ISOOppb. Methane concentrations and flight trajectory relative to source locations 
are illustrated in Figure [6} Methane concentrations in time are shown in Figure [7} Detail of plume 
extent for two intervals of flight trajectory (corresponding to the central region of Figure[6]) at 200m 
above ground level, and simulated concentration measurements at the aircraft, are given in Figure|8] 
Also shown are simulated concentrations as a function of relative northing for these trajectory 
intervals. Referring to pane (b), the longer right-hand tail of concentration with relative northing 
indicates that the aircraft's trajectory has a positive downwind component. Referring to pane 
(c), the source at approximately (20, 13)km contributes a concentration peak on the trajectory at 
relative northing of approximately 15km, and the source at approximately (20.5, 16)km a trajectory 
concentration peak at approximately 17.5km. The source at approximately (20.5, 15)km contributes 
a "shoulder" at approximately 16.5km. 

Inspection of Figures [6] and [7] shows that each source is upwind of at least part of the flight 
path, and that gas emanating from each source contributes to the concentrations. Each source 
is therefore identifiable in principle. Some sources (e.g. those near (20, 15)km) are close to a 
downwind section of the flight path, others (e.g. the source at (3,4)km) are relatively distant. We 
anticipate that the former will be more precisely modelled. Plumes from sources around (20, 15)km 
contribute to simulated concentrations on multiple downwind passes of the flight path, providing 
range information to resolve source location, whereas plumes from other sources (e.g. that near 
(20, 32)km) intersect the flight path just once. 

The survey area is partitioned into a 80 x 80 grid of 500m x 500m cells and a random field 
background model is assumed. Cell size is determined by a balance of physical (e.g. likely source 
extent) and computational considerations. The starting set of source locations and emission rates 
for RJMCMC is selected by sampling 15 locations from the initial gridded optimisation solution 
and weighting each cell by its estimated emission rate. RJMCMC is executed for 13000 iterations 
of which the first 3000 burn-in iterations are ignored, generating a dependent sample from the joint 
posterior distribution of parameters for source locations, half-widths, emission rates, background, 
measurement error and wind direction bias. Burn-in length was determined by inspection of trace 
plots and 10, 000 iterations post burn-in were judged to be sufficient to characterise the posterior. 



More formal procedures (e.g. Gelman-Rubin convergence diagnostic Gelman and Rubin (1992)) 
were not considered necessary. 

Estimated source emission rates are summarised in emission rate maps in Figure [9j Panel a) 
shows the initial optimisation solution, which identifies 9 of the 10 sources. Most estimated source 
locations are displaced downwind of the defined locations, closer to the flight path. The source at 
(3, 4)km, furthest from the flight path, was not identified. Panel b) shows the posterior median 
estimate. The MCMC median solution, in common with the optimisation solution, only identifies 
9 of 10 sources, though these are closer to their true locations. We summarise marginal spatial 
uncertainty in terms of the 2.5% and 97.5% credible values for source emission rates shown in 
Panels c) and d) respectively. 8 of the 10 known sources appear in Panel c). In Panel d) there 
are 3 spurious sources. For visualisation purposes only, source maps in Figure [9] a) , b) and c) are 
estimated from gridded source estimates generated at the end of each complete iteration (over all 
parameters) of the Markov chain. 



Figure 10 (a) compares background estimates from optimisation and MCMC (with credible in- 
tervals). Recalling that the actual background is constant at 1800ppb, both estimates are within 
O.lppb of the truth. Figure [lO|(b) compares unexplained residual concentration with simulated con- 
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Figure 6: Flightpath for the synthetic problem. Marker size and colour saturation represent strength and location 
of simulated methane concentrations. Arrow shows direction of air movement during flight. Locations and physical 
extent of defined sources (each with emission rate of 0.1 m'^s"^) are shown as crosses with overlaid black circles whose 
radii represent physical extent of the sources. 



centration from optimisation (red) and MCMC (black). For a good model fit, we expect residuals 
to be zero-mean and show no relationship to the measured concentration. The MCMC residuals 
are relatively well distributed around zero. For optimisation, residuals corresponding to simulated 
concentrations close to true background (ISOOppb) are small; for larger simulated concentrations, 
residuals are large and positive since the Laplace prior over source emission rate ^ penalises source 
strength, generally resulting in positive residuals. Source locations are constrained to the centres 
of grid cells for the optimisation solution, but not for the mixture model estimate. The interested 
reader should note that the case presented is a typical example from a number of simulated cases 
considered but omitted for brevity. 
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Figure 7: Simulated methane concentrations along flight path as a function of time for the synthetic problem. 



4.2. Landfills 

The analogous analysis procedure is adopted for modelling landfill measurements. For initial op- 
timisation, the survey area is partitioned into a 100 x 100 grid of 400m x 400m cells. The starting 
point for RJMCMC is chosen by sampling 5 locations from the optimisation solution, weighting 
each cell by its estimated emission rate. 



The estimated source emission rates are shown in Figure 11 Panel (a) shows the initial optimisation 
solution; sources within the landfill boundaries are supplemented by "ghost" sources downwind of 
the landfills. These are likely due to errors in the wind field data or inadequacies in the plume model 
as well as changes in the wind direction during the prolonged gas transit times to the measurement 
locations which are up to 15km away which cannot be incorporated within the plume representation. 
Panel (b) shows the posterior MCMC median source estimate. Sources are centred within each 
landfill with just a single additional "ghost" source downwind of the eastern landfill (28, 19)km, 
implying that correcting for wind direction bias (of approximately 2°) improves inversion. Panels 
(c) and (d) show similar characteristic to Panel (b). Interestingly, no other spurious sources appear, 
suggesting strong spatial localisation of sources in this case. To our knowledge there is no actual 
gas emission at the "ghost" source location in Panel (2) (which also appears in Panel (1)). If indeed 
the "ghost" source is an artefact, more sophisticated wind field corrections would be necessary to 
eradicate it. However, land fills are uncharacteristically strong emission sources that are detectable 
at extreme ranges with correspondingly extended gas transit times. For more representative sources 
transit times are much shorter. Figure [l2| (a) shows estimated background concentration along the 



flight path for initial optimisation and MCMC. Figure 12 (b) gives fitted residuals against measured 
concentrations for initial optimisation and MCMC. Residuals for higher concentrations are generally 
larger and positive, indicating underestimation of source emission rate. 



4.3. Flare stack 

The analysis procedure is similar to that adopted for the applications above. The survey area is 
partitioned into an 50 x 50 grid of 300m x 300m cells. The RJMCMC starting point is chosen 
by sampling 5 locations from the initial optimisation solution, again weighted by emission rate. 
Figure [3] shows a clear discrepancy between plume direction and mean wind direction predicted 
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Figure 8: Plume concentrations for the synthetic problem: (a) Detail of central region of figure [6] showing simulated 
concentrations along flightpath as black dots, size proportional to concentration. Blue dots are also used to indicate 
source locations. Background concentration is l.Sppm. Mean direction of air movement shown by the blue arrow, 
(b) Simulated concentration as a function of relative northing for left-hand pass, (c) Simulated concentration as a 
function of relative northing for right-hand pass. 



by UKMO wind field data. The mixture model, incorporating a constant wind direction bias 
parameter, successfully corrects this. Inspection of Figure [3] suggests a prior wind direction bias 
of approximately —18°. The corresponding posterior 95% credible interval is estimated to be 
[—18.12,-17.2]°. This uncharacteristically large wind bias is attributed to the flight being in the 
late afternoon (as the ABL subsides) and is situated at the coast where winds are inherently less 
predictable. 

Source emission rates in Figure [13] are estimated using corrected wind directions, otherwise the 
initial optimisation solution (Panel (a)) would be severely compromised. The posterior median 
mixture model result (Panel (b)) is very similar. Panel (c) shows that 2.5% credible values from 
the mixture model are < 0.004m'^s~^. Marginal 97.5% credible values in Panel (d) suggest greater 
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Figure 9: Source emission rate maps for the synthetic problem: (a) Estimate from initial optimisation, (b) Median 
estimate from mixture model, (c) Marginal 2.5% credible value from mixture model, and (d) Marginal 97.5% credible 
value from mixture model. For ease of comparison the mixture model results are presented on the same grid cell size 
as the optimsation solution. The locations of defined sources (each with emission rate of O.lm^s"^) are shown as 
black crosses. The locations and emission rates of estimated sources are colour-coded according to the scale. The 
synthetic flight path is shown as a red line. 



uncertainty in flare stack location (e.g. compared with the landfill), despite good initial optimisation 



and posterior median estimates. Findings from Figure [T4| are similar to those from Figure 12 



5. Discussion 



Detection and location of sources of gas emissions into the atmosphere is a topic of intense interest. 
In this work we describe a method for detecting, locating and quantifying such sources using 
remotely obtained gas concentration data. The method developed is broadly applicable to any gas or 
passively transported, detectable atmospheric component, such as aerosols, radon, smoke, ash, dust 
and viruses. The method can be applied using ground based or airborne concentration data collected 
using point or line-of-sight sensors. Here, atmospheric point concentration measurements are 
modelled as the sum of a spatially and temporally smooth atmospheric background, augmented by 
concentrations arising from local sources. We model source emission rates with a Gaussian mixture 
model and use a Markov random field to represent the atmospheric background concentration 
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Figure 10: Diagnostics for the synthetic problem: (a) Estimated background concentrations along the flight path as 
a function of time. The initial optimisation is shown in red and the mixture model is shown in black; median (solid), 
2.5% and 97.5% credible values (both dashed) shown, (b) Residual versus simulated methane concentrations from 
initial optimisation (red) and median from mixture model (black). 



component of the measurements. A Gaussian plume atmospheric eddy dispersion model represents 
the gas dispersion process between sources and measurement locations. An initial point estimate 
of background concentration and source emission rates is obtained using mixed ^2 — ^1 optimisation 
over a discretised grid of potential source locations. Subsequent reversible jump Markov chain 
Monte Carlo inference provides estimated values and uncertainties for the number, emission rates, 
locations and areas of sources, as well as atmospheric background concentrations and other model 
parameters. We investigate the performance of the approach first for a synthetic inversion problem. 
We then apply the model to locating and detecting of gas emissions using actual airborne data (a) 
in the vicinity of two landfills, and (b) in the vicinity of a gas flare stack. All analysis was performed 



using MATLAB (2011). 



As discussed in the introduction, individual model components and inference tools used in this 
work are commonplace in the modelling literature. The combination of components and tools, 
pulling together physical constraints with rigourous analysis, has proved useful for the remote 
sensing applications considered in this work. Nevertheless, to the best of our knowledge, this is 
one of the first applications of Bayesian inference using reversible-jump MCMC to simultaneous 
multiple source and smooth spatio-temporal background estimation. The Gaussian plume model is 
a particularly simple steady-state approximation to dispersion of a gas release into the atmosphere. 



widely used in the environmental modelling literature (e.g. Gifford (1976)). In this work, the 
plume model provides a reasonable basis for estimating known flare stack and landfill locations. 
However, we find that correcting bias in (predicted) wind directions supplied by UKMO improves 
inference. There is evidence, in the form of a "ghost" source downwind of the eastern landfill, that 
a simple bias correction is not adequate, and that a more sophisticated approach (e.g. a slowly 
varying wind-direction bias) might be beneficial. There is considerable opportunity to achieve this 
within the Bayesian modelling framework. Incorporating plume model uncertainty in the initial 
optimisation can be achieved in some sense by considering optimisation over a representative set 
of forward model matrices A (in ([l|), rather than a single choice. Predicated on the availability 
of wind field data of adequate quality we might also consider more sophisticated plume models, 
e.g. plumes following wind flow lines, or from computational fluid dynamics. In the MCMC case. 
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Figure 11: The estimated source emission rate maps for the landfill application: (a) Estimate from initial optimisation, 
(b) Median estimate from the mixture model, (c) Marginal 2.5% credible value from the mixture model, and (d) 
Marginal 97.5% credible value from the mixture model. All dimensions in km. Emission rates in m^s~^. Each panel 
shows a common subregion of the original 40fcm x 40fcm domain (referenced with respect to the origin) within which 
all sources are estimated. For ease of comparison the mixture model results are presented on the same grid cell size 
as the optimsation solution. Polygons indicate the perimeters of the landfills. The fiight path is shown as a red line. 

we assume a-priori that sources can occur with equal probabiUty at any location. For RJMCMC 
we sample those grid locations from the optimisation solution with the greatest emission rates 
as a starting solution. Work continues to explore incorporation of spatial prior distributions for 
source location, for example using Polya trees to encode some degree of source clustering. There 
is scope to develop more sophisticated background models incorporating parameters known to 
influence background methane concentration (such as topography). It would be interesting to 
explore modelling background as the superposition of a number of distant sources. Our field 
experience suggests that natural gas seeps can be intermittent, requiring adaptation of our model 
formulation. Smoothly varying gas release rates could be accommodated relatively simply. 

The Gaussian plume model provides an elementary means of modelling gas transport from source to 
measurement location under ideal steady state wind field assumptions, allowing rapid estimation of 
forward model matrices A at the expense of accuracy and precision. Given inherent uncertainties 
in the estimates of wind field parameter values supplied by UKMO, we consider the Gaussian 
plume adequate for the purposes of the current work. For example in the landfills application, 
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Figure 12: Diagnostics for the landfill application: (a) Estimated background concentrations along the flight path as 
a function of time. The initial optimisation is shown in red and the mixture model result is shown in black; median 
(solid), 2.5% and 97.5% credible values (both dashed) shown, (b) Residual versus measured methane concentrations 
from the initial optimisation (red) and median from the mixture model (black). 



assuming ideal wind field conditions over an interval of approximately 10 minutes for wind speed of 
approximately 6.5ms~^, suggests that the plume model is appropriate for measurement locations 
within approximately 4km downwind of a source. More distant measurement locations require ideal 
conditions over longer periods for the Gaussian plume model to be appropriate. Nevertheless, even 
over larger distances, the Gaussian plume is likely to provide reasonable approximation to reality 
provided the wind field remains relatively steady. 

As implemented in the current work, optimisation is used to provide an initial point solution for 
inversion on a spatial grid. Subsequent Bayesian inference gives a more flexible grid-free mix- 
ture model framework within which to estimate the joint posterior distribution of all parameters, 
providing in particular estimates for parameter uncertainty. Early attempts at inversion followed 
a stepwise approach in which atmospheric background was estimated prior to, and independent of 
emission sources. The current approach, involving simultaneous estimation of background, sources 
and wind fleld characteristics improves performance. We also explored Bayesian inference on the 
same spatial grid used for the initial optimisation. The very large number of potential source 
locations makes this computationally intensive. 

Our experience of processing multiple survey datasets has made clear the need for rigorous data 
management and pre-processing procedures, e.g. in the merging of spatio-temporal data from 
different sources (e.g. aircraft and UKMO). Efficiency of inference can be improved for a given 
deployment by specifying a flight trajectory (or sequence of flight trajectories) appropriately, given 
prevailing wind conditions and prior information concerning likely source locations. Methods of 
statistical experimental design are central to achieving this for both airborne and ground based 
line-of"Sight gas sensors. It is also strongly desirable to have a means of confirming the quality 
of inference, particularly of source location and release rate, using a persistent known gas source 
within the region of interest. In some cases, this might take the form of an existing methane source 
(such as a flare stack or landflll), or perhaps a small controlled release (e.g from a gas cylinder). 
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Figure 13: Estimated source emission rate maps for the flare stack appflcation: (a) Estimate from the initial optim- 
isation, (b) Median estimate from the mixture model, (c) Marginal 2.5% credible value from the mixture model, and 
(d) Marginal 97.5% credible value from the mixture model. All dimensions in km. Emission rates in m'^s^^. Each 
panel shows a common subregion of the original 15fcm x 15fcm domain (referenced with respect to the origin) within 
which all sources are estimated. For ease of comparison the mixture model results are presented on the same grid 
cell size as the optimsation solution. The black annulus indicates the location of the flare stack, a point source. The 
flight path is shown as a red line. 
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Figure 14: Diagnostics for the flare stack application: (a) Estimated background concentrations along the flight path 
as a function of time. The initial optimisation is shown in red and the mixture model result is shown in black; median 
(solid), 2.5% and 97.5% credible values (both dashed) shown, (b) Residuals versus measured methane concentrations 
from the initial optimisation (red) and median from the mixture model (black). 



Appendices 

A. Background model 
A.l. MRF background model 

We model the background as a Gauss-Markov random field, a class of graphical model. Markov 
random fields are joint distributions for variables Xi, . . . ,X]sf and an associated undirected graph 
G = (£, V). The vertex set V := {1, . . . , A^}, and the edge set £ is a subset of V x "V. The graph 
specifies the conditional independence structure of random variables as follows. If three sets of 
variables A,B,C C V are such that the set B separates A from C in the graph G, then (Xj)jg^ 
must be conditionally independent of {Xi)ii^c given (Xi),^^. A simple special case of a Markov 
random field is a Markov chain, wherein G is simply a linear graph. In a Gauss-Markov random 



field, the random variables are also assumed to be jointly Gaussian. It can be shown (e.g. Speed and 



Kiiveri (1986)) that for Gaussian variables, the conditional independence condition is equivalent to 
the precision matrix J being sparse with respect to the graph G. That is, for i different from j, 
Jij 7^ if and only if there is an edge between node i and node j in G. 

For the random field model, we particularise (|4]) as: 

/(/3)ocexp{-i/3%/3} (12) 



where again is sparse with respect to a graph to be designed. As stated in Sec. 3.1, we take P 
to be the identity matrix, so that our background estimate b is simply equal to (3. 

In designing a Gauss-Markov random field to model the background field, we seek to capture 
two effects. First, the background should be smooth. Second, since the background concentration 
travels with wind, it should be smoother along the direction of the wind. To model these two effects, 
we introduce two different types of edges in the graphical model. The first are edges connecting 
adjacent measurements vertices. These ensure overall smoothness of the background with respect 
to time and space. The second type of edge concerns the wind. At each measurement point, we 
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consider the line along the wind direction from this point. We find the next measurement point 
that crosses this line, and connect an edge in the graph between these two points. Thus the second 
point is as close as possible to directly in line with the wind from the initial point. "Wind-linked" 



points along the trajectory for the landfill application is given in Figure 15 
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Figure 15: Illustration of "wind-linked" points along the flight trajectory for the landfill application 



The graphical structure outlined above determines the sparsity pattern of J^, but not the values 
themselves. The structure of Jp is as follows: 

(ij)e£ 

where Ajj is a matrix non-zero only at the intersections of the ith and jth rows and columns, and 
the non-zero elements of Ajj are given by: 

1 -1 " 
-1 1 

This may be alternately stated as, for any vector x: 

(ij)6£ 



X 
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The values aij represent the strengths of the hnks, and they are chosen as fohows. We assume 
that if the background value were measured in the same packet of air at two different times, the 
change in background concentration would roughly grow with the time difference. Similarly, at 
two different points in space at the same time, the difference in background concentration would 
grow with the spatial separation. This suggests that the strength of each link should be determined 
by the time difference and the spatial distance between the two connected points. In particular, 
we take as parameters two constants ct and cd- The parameter ct is the expected change in 
background concentration per second difference in measurement time, and C£) is the change in 
background concentration per meter separation. The strength Oij is then given by: 

1 

an 



{ct at + CD Aoy 

where AT is the difference in time between the two measurements and AD is the spatial distance 
between the two measurements. For the links between adjacent measurements, we calculate AD 
as simply the distance from first measurement point to the second measurement point. For the 
wind-based cross-links, we use the distance from the second measurement to the point in space 
where the packet of air at the first measurement would be, assuming that it travels with the wind 
at the wind speed and direction as measured at the first point. Other parametric forms for aij, 



e.g. motivated by diffusion-advection transport (see A. 2 below), could be considered. 



A. 2. Polynomial background model 

Based on the idea that the background field is slowly varying in space and time, we can also 
construct a basis for a smooth, low-dimensional manifold using low-degree Chebyshev polynomials. 
If desired, it can be trivially extended to more dimensions to account for regression on other 
variables, such as pressure, temperature, etc. 

The background field is modelled as a function of space x, y, z and time t, and coefficient vector /3: 

b{x,y,z,t; P) = ^ l3ij^k,iTi{x)Tj{y)Tk{z)Ti{t) 

i,j,k,l 

where Tn is the Chebyshev polynomial of the second kind of order n, with the coordinates properly 
transformed such that each dimension x,y,z, and t are mapped onto the interval [—1, 1]. 



The unsealed Chebyshev polynomials of the second kind Arfken ( 1985 ) are defined by the recurrence 
relation (x £ [—1, 1]) 

Uoix) = 1, Ui{x) = 2x, Un+lix) = 2xUn{x) - Un-l{x) . 

One of the advantages of the Chebyshev polynomial model is that it has an analytical form, and 
therefore it is very easy to evaluate the function and its derivatives on any point in space, particu- 
larly those outside the flight track. 

The regularisation and conditioning of this type of background consists of three components: 

Regularisation with respect to a reference background. A constant reference background 
level bo = PPq can be used for regularisation. Penalisation for deviation of b from bo can be 
written as the quadratic form 

i?i(/3) = (/3-/3of/(/3-/3o). 
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Smoothness regularisation. We choose as a smoothness measure the sum of squares of the 
curvature of the background field evaluated at G = {^j = {xi,yi, Zi,ti)}i, a discretisation of 
the spatio-temporal domain. 



R2{P) = E 



dH\^ , /a26\2 /a26\2 /525\2 



+ 



+ 



(13) 



and the analytical form (A. 2) allows us to easily write (13) as a quadratic form: 

= (/3-/3o)^J2(/3-/3o) 

where J2 is a symmetric semi-definite matrix. 

Wind transport regularisation. The background field can be made to obey the transport equa- 
tion along the wind field w(x, y, z, t\. 



dh 
Wt 







Similar to (13), we formulate the wind transport penalty as a quadratic form 

2 



dh 

w(Ci)-v,,2,,,6(ei) + — (e^) 



RM = E 

6 eg 

= (/3-/3o)^J3(/3-/3o), 



where J3 is a symmetric semi-definite matrix. 
In the end, the three penalisations can be added as a single quadratic form 

(P - /3o)^J(/3 - /3o) ) where J = ml + /i2^2 + fJ-iJs , 

and ^1, /i2, and 113 are relative weights that allow an expert to choose the strength of each regu- 
larisation term. 



B. Initial parameter estimation 



We solve the optimisation problem ^ based on an alternating direction method. The objective 
function is not entirely separable with respect to the variables s and /3, since they are coupled 
through the data fitting term \\As + P(3 — y|p, but the rest of the terms and the constraints are 



uncoupled. These type of split methods were devised originally (Pukushima (1992)) for large-scale 



separable problems, but variants have gained a lot of attention recently to solve other problems, 
especially in the fields of image processing, machine learning and compressed sensing (Yin et al. 



(2008); Yang and Zhang (2009); Deng and Zhang (2011)) where the objective function can be 
written as the sum of convex subproblems, each with special structure and characteristics that can 
be exploited, such as the sum of smooth and non-smooth terms (like the £2-^1 problem that we 
have). Even if the optimisation problem can be written as a standard linear or quadratic program 
-and their respective theories and solvers are very mature- these new methods can have better 
practical and computational properties. 

In our case, the main optimisation problem ([s]) is indeed a quadratic objective function with linear 
inequality constraints, and J is a positive semi-definite matrix, so it could be solved with standard 
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convex quadratic solvers. However, the problem has a special structure that can be exploited and 
we have found that a iterative split technique converges much faster than using Matlab's quadprog 
routin^B 

The algorithm can be described as an alternating sequence of optimisation of the background and 
ground sources. In the first stage, ([s]) is minimised only with respect to (3, the parameterised 
background, while keeping the current estimate of s fixed. Analogously, in the second stage, ([s]) 
is minimised only with respect to the ground sources emission rate s, while keeping the current 
estimate of the background fixed. Then, a criterion of optimality is computed and, if satisfied, the 
current solution is deemed to be close enough to the global minimum. If it is not satisfied, the 
algorithm iterates again through the first and second stages. 



Step 1: Estimate the background. Although it is physically reasonable, the non-negativity 
constraint of the background, < P/3, is actually not enforced explicitly, since in practice it 
is never active. The sub-problem 



min 2^ II As + P(3- yf + f (/3 - l3ofj{(3 - (3^) 

f3,-W e 



subject to Pp < y + T . 



is solved by the method of augmented Lagrangian for inequalities (Nocedal and Wright 



(2000)). To do this, the inequality constraint is converted into an equality constraint by 



introducing a slack non- negative variable w: 



min ^\\As + Pp- y||2 + ^{(3 - PofJi(3 - (3, 



subject to 



-PP + y + r + w = 
w > 



This, in turn, is solved by moving the equality constraint into the objective function and 
introducing its corresponding Lagrange multiplier z G M", to obtain 

n 

min L(/3, z;ri) = ^ \\As + P(3 - yf + f (/3 - Pofj{l3 - (3q) + Y^^HP)^ 
subject to w > 
where Cj is the i-ih. component of the constraint 

c = -P(3 + y + T + w 
and tp is the auxiliary function defined as 

1—^6 otherwise. 

and this problem is solved with a combination of Newton's method and projected gradient 
to maintain the feasibility w > 0. For more details on the practical considerations with this 



type of methods, see §17.4 of (Nocedal and Wright (2000)) 



^Matlab 2008a 
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Step 2: Estimate the sources. The subproblem 

mill tAII^s + P/3 - ylP + AHQslli 

subject to < S < Smax 

has simple bound constraints and we use a gradient-projection method to solve it. We use a 



majorise-minimise (Chouzenoux et al. (2011)) method to solve for this quadratic form with 



simple bound constraints < s < s„ 



min ^\\As + PI3- + ^(/3 - /3o) Jil^ " /^o) + W s 



Note that the term ||(5s||i has been replaced by q^s, where q = diag((5). This is due to 
the fact that since we have that ||s||i = l"^s when s > 0, and therefore the term is fully 
differentiable in the half-space s > 0, and no concerns of the non-differentiability of the 



•norm must be taken into account. 



C. Mixture model 
C.l. MCMC step types 



Referring to section 3.3, we write the parameter set as 6 for brevity. 6 consists of source parameters 
z, w, s, background parameters /3, measurement error standard deviation and (potentially) wind 
direction and other plume bias and uncertainty terms. 

Metropolis Hastings 

In conventional MCMC using the Metropolis Hastings algorithm, we construct a Markov chain with 
stationary distribution corresponding to the posterior distribution f{6\'D) given observed data V 
consisting of observed airborne concentrations y (and potentially wind field measurements). For 
transitions between two states and {0^,0^}, we ensure reversibility of the Markov chain 

by imposing the detailed balance condition: 



/ 



f{9^\v,eT,)q{e^,e',\e^)a{e^,e',\e^)de^ = / f{e'jv,eT,)q{e',,e,\e^)a(ele^\ej,)de', 



for acceptance probability a and proposal distribution q, from which (11) emerges. In the current 
work, the proposal distribution q always corresponds to a random walk such that (jf(0K, ^kI^k) = 
q{6'i^,9i^\9-f=:) (e.g. we might use 0^ = 0^ + e for (multivariate) Gaussian e). The manner in which 
the chain explores the posterior distribution for dependent variables or multi-modal distributions 
(such as our mixture model) can be improved by updating subsets of parameters in a particular 
order and adopting a good starting solution. 

Reversible jump Metropolis Hastings 

In RJMCMC, we construct a Markov chain which satisfies an extended balance equation: 

J f{o^\v,eT,)g{(t}\eT,)j{e^\e^)a{e^,e'jeT,)do^d^ = J f{o',\v,eT,)g'{ct)'\e^)j{e',\e^)a{o',,e^\0T,)de'j(t,' 
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where and 0^ are now of different dimensions . We facilitate dimension-jumping by constructing 
augmented sets of variables {6^, 4>} and {9'^, cf)'} which are of the same dimension, using specified 
bijective functions h to move between them, such that (0'^, cf)') = <P) (and (0^, cf)) = h'{0'i^, cf>'), 
with h' = ). We also specify the joint distribution g{cf>\0-f^) (typically uniform or Gaussian; 
g'{cf>'\6j^) is known given g and h). j{6fi\6^) is the probability that a particular dimension- jumping 
move will be attempted from state {0^, ^k}- 

The Metropolis-Hastings acceptance probability a(0K, ^'^I^k) (corresponding to ([XT])) for a dimension- 
jumping move becomes: 



a(6>«,6>'^|6>«) =maxh 



fie'jv,e^)jie',\e^)g'icf)'\e^) 
f{o^\v,e^)j{o,\o^)g{cf)\e^) 



d{e'^,cf>') 



d{9^,cf)) 



(16) 



where the Jacobian is easily evaluated given h. The acceptance probability for the reverse move is 
similarly: 



a{0'^,0^\e^) = max\ 1, 



f{e^\v,o^)j{o^\o^)g{cf)\e^) 
f{e'\v,en)j{e'^e^)g'{ct)'\e^) 



d{e,cf)) 



d{e',cf)' 



We use RJMCMC to update the number of sources m by independent birth-death and coalesce-split 
steps. 



An illustrative split step 

We include a brief description of a typical split step to illustrate RJMCMC. In the split step, we 
select an existing source j* at random and consider splitting it in two. In this case, 0^ is the triplet 
{zj*, Wj*, Sj*}. We draw random variables {r^,r^,rs} from uniform distributions on appropriate 
intervals with zero mean, which to create two new sources j*"*" and j*^ to replace j* . The new 
source parameters are given by: 



z,.± = z,.±r„ r.(l)~[/([-^,^]), r.(2)~t/([-^,^]) 

Wj^± = Wj^±ry,, ryj U{[-Eyj/2,Euj/2]) 
Sj*± = Sj* ±rs, rs ^ U{[-Es/2,Es/2]) 



In this case, the Jacobian of the transformation (from (|16|)) is equal to 2 x 2 x 2 = 8. If we assume 
that source location components, widths and emission rates have independent uniform priors on 
[0,-R2i], [0,-^22], [0,Rw] and [0,Rs], expressions for f {6'i^\T> , 6,^) and f {6 k.\T> , 6,^) emerge: 

f{e^\V,e-^) ~ f{V\e^,e^) x p ^ p p x /(Prior for all other parameters) 



f{e',\v,e^)r. f{v\e',,e^)x 



1 



RzlRz2RwRs 



X /(Prior for all other parameters) 



Noting that g = {EziEz2EwEs) ^ and that g' = 1, using (16) we obtain the acceptance probability 
for splitting source j*: 
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a (split source i*!^^) = max < 1 



/(^kI^)^k) jX^kI^k) RziRz2RwRs 



X 8 



} 



in which jump probabilities j{6\0K) and j(0'|07t) are equal. 
C.2. MCMC procedure 

A new state of the Markov chain having the posterior distribution of model parameters given 
observed data as stationary distribution, is created from the current state of the chain using the 
following procedure: 

Update source parameters, {z,w, s}. Firstly, source locations are updated sequentially using 
a Gaussian random-walk. Candidate locations outside the spatial domain of interest are 

rejected. Secondly, source widths are updated sequentially using a Gaussian random-walk. 
Finally, source emission rates are updated sequentially using a Gaussian random-walk. 

Update measurement error, cry. Measurement error variance is updated using a Gaussian ran- 
dom walk on logg ay- 

Update background parameters, /3. The background parameter vector /3 is updated by sampling 
from its full conditional. 

Optionally update wind direction bias, and plume opening angles, ujh, i^y- The bias in 
wind direction is updated using a Gaussian random walk modulo 27r. Then the plume opening 
angles arc updated if required using a Gaussian random walk on the log scale. 

Update number of sources, m. With equal probabilities of 0.25, one of 4 candidate states is 
generated, which if accepted, either increases (by birth or splitting) or decreases (by coales- 
cence or death) the number of sources by 1. Steps involving new states in which m < or 
m > mmax are rejected, where mmax is an upper bound for the number of sources. 
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